Monitoring spring phenology of individual tree crowns using drone‐acquired NDVI data
Abstract
Quantifying the timing of vegetation phenology is critical for monitoring and modelling ecosystem responses to environmental change. Phenological processes have been studied from landscape to global scales using Earth observing satellite data, and at local scale by in situ surveys of individual plants. Now, data acquired from multi‐spectral sensors on drone platforms provide flexible opportunities for monitoring phenology from individual plants to small ecosystem scales efficiently, allowing community and species level information to be derived. We captured a time‐series of drone‐acquired normalized difference vegetation index (NDVI) data with a multi‐spectral sensor (Parrot Sequoia, (Parrot, France)) over a highly heterogeneous ecosystem in Cornwall, UK, during a period of spring green‐up. We monitored NDVI trajectories at the individual crown and species’ level. For deciduous crowns, we derived metrics representative of spring phenological stages: Start‐of‐spring (SOS), middle‐of‐spring green‐up (MOG) and start‐of‐peak greenness (SOP) using a logistic function. While the exact timing of SOS, MOG and SOP appeared susceptible to understorey effects and saturation of the NDVI, relative timing of green‐up for a subset of species was plausible in relation to phenological observations from an extended geographic region and in situ plant area index (PAI) measurements. In evergreen vegetation (Pinus spp.) subtle changes were also detected through the growing season. The impact of illumination differences was analysed for image pairs during leaf‐off and leaf‐on conditions. While significant, these effects were small (mean absolute NDVI deviation of up to 0.034 for leaf‐off, 0.013 for leaf‐on conditions), meaning that data captured under both constant direct and diffuse irradiance conditions can be used together and that cloudy conditions should not lead to data gaps. We conclude that the capability of drone‐mounted multi‐spectral instruments for spatio‐temporal characterization of crown‐level phenology shows great promise for improving the understanding of intra‐ and inter‐species differences in strategy, and offers an efficient means of doing so over areas of a few hectares.
Introduction
The advantage of centimetre scale spatial resolution remote sensing in the vegetation context lies within its ability to resolve the fundamental ecological unit: individual plants (Marconi et al, 2019). In mixed species ecosystems, resolving individual plant crowns is necessary if seeking to analyse species‐specific phenological trends, plasticity and responses to extreme events, for example the earlier green‐up of some species in response to shifts in temperature and consequently growing degree days (Vitasse et al., 2010), or the expression of stress (Adams et al., 2015). It has previously been concluded that satellite‐based data are not suitable for monitoring the phenology of individuals or species (Polgar & Primack, 2011). A wealth of satellite‐based studies on tree phenology has predominantly focused on the ecosystem, landscape or forest stand scale (Melaas et al., 2013; Pastor‐Guzman et al., 2018; Pennec et al., 2011; Walker et al., 2012). While studying species‐level phenology using satellite data is possible, it requires larger mono‐species regions due to spatial resolutions which are too coarse (>10 m) to resolve individuals in mixed species systems (Berra et al., 2019; Delbart et al., 2005; Han et al., 2013; Liu et al., 2017). In order to investigate species‐specific phenology, to relate this to local conditions as well as quantify the effect of landscape heterogeneity on large‐scale satellite‐derived metrics efficient monitoring at the individual level is required (Klosterman et al., 2014). The lack of ecosystem‐scale ground phenological data which can serve as reference over adequately large spatial extents (e.g. 250 m pixels) has previously been highlighted as a major limiting factor for the validation of satellite‐derived estimates (Cleland et al., 2007). The use of fine spatial resolution satellite data (e.g. WorldView, PlanetScope) is limited due to their financial cost and missing observations in areas with frequent cloud cover (Bolton et al., 2020; Hufkens et al., 2012). The remote observations which best enable phenological monitoring at the species and individual level are cameras mounted in situ, such as the global network of phenocams providing daily measurements (Brown et al., 2016; Polgar & Primack, 2011). However, phenocam images suffer limitations due to oblique viewing angles such as measurement saturation at low leaf‐area‐index (Keenan et al., 2014) and canopy occlusion resulting in limited spatial coverage.
Recently, the acquisition of crown time‐series data over areas of multiple hectares in extent has become a possibility through unmanned aerial vehicle (UAV or drone) platforms. Drones provide advantages over phenocams by covering larger areas with close to nadir view angles enabling calibrated reflectance‐based metrics (Piao et al., 2019) and unlike optical satellite instruments can acquire useful data during overcast conditions. Coverage limitations of drones (approx. 5 ha according to a cost‐benefit analysis of acquisitions by multi‐rotor (Matese et al., 2015)) mean that drone‐based surveys are best suited to capture data over plot scales that are broadly representative of a larger region, for example exhibiting similar species composition and climatic conditions. As compared to surveys using piloted aircraft or commercial satellite imagery, drone‐based acquisitions can be lower cost (dependant on the size of the study site, accessibility and cost of labour) with the only limitations to revisit times for data acquisitions being meteorological conditions and accessibility of the study site (Duffy et al., 2017).
Drone‐based monitoring of phenological transitions of individual tree and shrub crowns has recently been proven possible (Berra et al., 2019; D’Odorico et al., 2020; Klosterman et al., 2018; Park et al., 2019). Most such studies use primarily RGB (red, green and blue band) photography from consumer‐grade cameras (Berra et al., 2019; Klosterman et al., 2018; Park et al., 2019). Supplementing RGB with information from the near‐infrared (NIR), where healthy green vegetation is highly reflective allows the derivation of plant indicators such as the normalized difference vegetation index (NDVI) which shows the presence and measures the photosynthetic capacity of the canopy (Tucker, 1979). It has previously been shown that NDVI remains more sensitive to leaf development than RGB‐based indices (Brown et al., 2017). NDVI from modified consumer‐grade cameras proved poor at crown scale (Berra et al., 2019), but compact radiometric, multi‐spectral sensors are now widely available for lightweight drones and allow for simultaneous narrow‐band NIR measurements enabling potentially, more robust NDVI‐based monitoring of changes in vegetation (Aasen et al., 2018; Duan et al., 2017), and which can also be customized to monitor evergreen species phenology using pigment indices (D’Odorico et al., 2020).
Data from drone‐based multi‐spectral sensors thus appear promising for the reliable monitoring of leaf development from budburst to senescence and the subsequent derivation of phenological stages. Knowledge of the robustness of this data and its sensitivity to changes over time is however crucial and as yet, untested. Specifically, the geometric consistency (Marques et al., 2019) and the reproducibility of vegetation products under varying illumination conditions (Damm et al., 2015; Gamon et al., 2006; Lee & Kaufman, 1986) should be assessed, which leads to specific questions that are posed in this study.
We investigated the robustness of the NDVI time‐series and sought to extract phenological indicators of spring green‐up at the individual level from drone‐based radiometric, multi‐spectral data captured through a spring green‐up time‐series. In particular, we set out to answer the following research questions:
- Are NDVI orthomosaics derived from multi‐spectral drone imagery sufficiently robust over time to track changes at individual crown level?
- Is the georeferencing accuracy following standard georeferencing procedures sufficiently high to ensure data are co‐registered well enough for time‐series reconstruction for individual crowns?
- How large are NDVI variations due to different atmospheric conditions during acquisition?
- Can the drone‐derived NDVI data reveal phenological differences at the species and individual crown level?
- What species‐specific features can be identified from the NDVI time‐series?
- How plausible are derived phenological metrics (start‐of‐spring, middle‐of‐spring green‐up and start‐of‐peak greenness) and what do they reveal about intra and inter‐species variability?
Answering these questions represents an important step towards the operational drone‐based monitoring of phenology at the individual crown level.
Materials and Methods
Study site
The study site ‘Trelusback Farm’ is located in West Cornwall, UK (50°12'10.0"N 5°12'22.6"W) with the studied area covering 4.5 ha. It is managed for the purpose of wildlife conservation and includes a small mixed woodland of predominantly sweet chestnut (Castanea sativa), sycamore (Acer pseudoplatanus) and alder (Alnus glutinosa) as well as isolated free‐standing trees and shrubs, mainly hawthorn (Crataegus monogyna), oak (Quercus spp.), grey willow (Salix cinerea) and gorse (Ulex europeaus). The site was selected due to its tree species richness, accessibility and due to being a heterogeneous system with variation in canopy density which is a common feature of many natural and semi‐natural ecosystems worldwide.
Drone data acquisition and processing
Data were acquired between 22 March 2019 and 23 July 2019 using the Parrot Sequoia multi‐camera array (Parrot, France), which possesses four cameras with different band‐pass filters in order to record light in the green (550 nm, 40 nm bandwidth), red (640 nm, 40 nm bandwidth), red‐edge (735 nm, 10 nm bandwidth) and NIR (790 nm, 40 nm bandwidth) wavelength regions. The Sequoia also records irradiance in equivalent bands with an upward‐facing irradiance sensor. The camera was mounted on a 3DR‐Solo lightweight quadcopter. For multi‐spectral data acquisition, the drone was flown in a pre‐programmed north‐south lawnmower pattern at 5 m/s ground speed and 70 m height in order to be able to cover the entire study site within a single flight with approx. 9‐minute acquisition time. This resulted in 8.1 cm mean ground sampling distance and >80% frontal and lateral image overlap as advocated by other studies using the Parrot Sequoia sensor (Assmann et al., 2019; Tu et al., 2018). Visual reference datasets for crown discrimination were acquired using a consumer‐grade RGB camera (Ricoh GRII), mounted on the same platform and flown at 6 m/s ground speed and 60 m height. Six ground control points (GCPs) were deployed in the field (see fig. 1 for reference), distributed following recommendations by James and Robson, (2014) to minimize systematic error within the photogrammetric processing. Two GCPs were redeployed in slightly different locations on 19 April 2019 and 19 May 2019, respectively, due to sub‐optimal coverage and plant over‐growth. Due to the limited number of GCPs, all were used within the photogrammetric processing workflow to guarantee the best possible georeferencing of the output. For an independent investigation of georeferencing accuracy, a subset of three datasets (22 March 2019, 19 April 2019, 23 July 2019) were assessed by excluding a GCP as check‐point and processed using the Pix4D (version 4.4.12, Pix4D, Switzerland) photogrammetric software. Pix4D was then used to generate a multi‐band orthomosaic for each date from single band images, using the ‘Ag Multispectral’ processing template. Noise filtering and medium surface smoothing were enabled. Keypoints image scale was set to full with a targeted number of 10'000 keypoints. Half image scale was used for densification and point density was set to ‘optimal’. Pix4D applies camera‐specific corrections (lens distortion and vignetting), the conversion of digital numbers to pseudo‐radiance (values homogenous to actual radiance (Parrot, 2017)) and the further calibration to surface reflectance using reference images of a levelled 50 x 50 cm panel of 44% nominal reflectance across all four bands of the Sequoia sensor (MosaicMill, Finland) taken prior to or post‐flight (see Fawcett and Anderson, (2019) for more details on the calibration method). The panel reflected pseudo‐radiance in the red and NIR bands was derived from calibration images and reported in Tab. SI1 as an indicator of irradiance conditions. Data from the irradiance sensor were used to eliminate possible small differences for data acquired during diffuse conditions but could not be applied to two diffuse acquisitions following 11 July 2019 due to malfunctioning of the irradiance sensor connection. The impact of this is expected to be small due to homogeneously overcast skies. The described workflow allows a more straight‐forward calibration of the data to surface reflectance compared to modified off‐the shelf digital cameras (Holman et al., 2019).
The NDVI as a metric of choice for vegetation monitoring over time is merited due to its sensitivity to the physical parameters of interest as well as greater consistency compared to common metrics derived from uncalibrated RGB photography (Yao et al., 2017).
Field reference data
A sample of individual trees (206 total) was manually digitized from the RGB orthomosaic and the species determined following a GPS‐assisted field survey and visual interpretation of the RGB data. A 1.5 x 1.5 m concrete surface within the study site served as unvegetated reference to control for erroneous NDVI variability over the duration of the study (Fig. 1).

To investigate the relationship of crown NDVI values with plant area index (PAI) during spring green‐up, a transect of digital hemispherical photographs (DHP) through the woodland and DHPs at 10 locations below individual crowns were taken on 10 dates from pre budburst of most species to full leafing out. Six dates corresponded directly to drone acquisitions during diffuse conditions and the maximum difference in time between drone and DHP acquisition dates was 5 days. DHPs were taken using a Nikon D7000 with a Sigma 4.5 mm F2.8 EX DC HSM Fisheye lens, along north‐south transects with 5 m spacing during periods of no direct sunlight (mostly overcast, some with scattered clouds). For individual crowns, DHPs were taken close to the stem and the image sector containing the stem was masked. Each image was processed using the Hemisfer software (version 2.2, WSL, Switzerland), applying the LAI‐2000 algorithm as implemented by Li‐Cor (USA), and the Li‐Cor recommended adjustments for individual crowns (Cutini & Varallo, 2006), detailed in the supplementary information. Areas of interest with the transects at their centre of approx. 17 m width and 75 m total length for the woodland and manually delineated individual crowns were extracted from the NDVI raster datasets to yield a time‐series of NDVI mean and standard deviation. For direct comparison with PAI observations of non‐matching dates, linear interpolation of the NDVI observations was performed.
Dates of budburst and emergence of the first leaf were extracted for 7 of the analysed species from the Woodland Trust's ‘Nature's Calendar’ which displays data collected by the UK Phenology Network (UKPN), for comparison with identified phenological metrics. This analysis was constrained to the counties of Cornwall and Devon (South West England) and the same year (2019) to account for climatic variations.
Assessment of crown NDVI values in direct versus diffuse conditions
At opposite ends of the time‐series, during leaf‐off and leaf‐on conditions, repeat acquisitions during contrasting illumination conditions with only 1‐2 days time difference were acquired (22 March and 24 March, 31 March and 1 April, 22 July and 23 July 2019). These were used to assess the impact of differing illumination conditions on retrieved crown NDVI.
NDVI density distributions for all crown pixels were generated. Furthermore, the coefficient of determination (R2), mean absolute deviation (MAD) and bias between the direct and diffuse NDVI values were compared for different crown aggregation metrics. The metrics assessed were the mean, median, and the mean after masking all but the 20% brightest pixels in the green band, per crown, a recommended method for mitigating the influence of the crown background (Berra et al., 2019). A two‐sample t‐test was performed with the null hypothesis that the means of the direct and diffuse per‐crown NDVI values were identical.
Green‐up modelling and phenophase estimation
Where through
are fitting parameters representing
: the value of t at the inflection point,
: the growth rate,
: the amplitude of increase and
: the lower asymptote. (eq. 2) was fitted to the data using the ‘nls’ function in the R stats package. After successful fitting of the model, the minima and maxima of the rate of change of curvature, calculated using the first and second derivatives of (eq. 2) could be used to identify SOS, MOG and SOP (Klosterman et al., 2014) (Fig. 2).

In a preliminary analysis, it was observed that the flowering of hawthorn resulted in a period of low NDVI within the time‐series which led to failure during model fitting. As the flowering event does not represent the spring green‐up modelled here (obscuring the signal of leaves), four hawthorn observations during this period (May to mid‐June) were excluded from the data used to fit the logistic function.
Results
Surface reflectance and NDVI orthomosaics
Twenty surface reflectance orthomosaics (Fig. 3, irradiance conditions for each are given in Table SI1) were generated with a mean reported georeferencing error of 0.023 m (standard deviation: 0.012 m). When using one GCP as check‐point, assessed for three datasets, the mean horizontal error (XY) was 0.023, whereas the mean vertical error (Z) was 0.149. As the horizontal error was lower than 1 GSD, the orthomosaic products had minimal visual spatial offset.

For the cloud‐free acquisitions pre green‐up, there were a lack of identified tie‐points within sparse tree crowns and an increased number of points identified on the ground below resulting in an incorrect representation of some crowns within the generated digital surface model (DSM) and consequently in the NDVI orthomosaics due to mosaicking issues. This could not be mitigated by disabling surface smoothing and noise filtering in Pix4D. For overcast acquisitions and post‐green‐up, crown shapes were found to be highly consistent.
NDVI sensitivity to illumination and aggregation metrics
Visual evaluations of surface reflectance products resulting from acquisitions during contrasting direct and diffuse illumination conditions showed that direct illumination resulted in cast shadows, either from the crowns on the ground or within crowns due to their structure. These shadowed areas were not as apparent in the NDVI product as it is formed by band ratios but fine‐scale variations were visible within single crowns if scaled appropriately (Fig. 4).

The density distributions of all crown NDVI values for two leaf‐off and one leaf‐on snapshot (Fig. 5) revealed that direct illumination conditions led to a larger spread in values both pre‐ and post‐green‐up. The differences for leaf‐off conditions were assessed using two pairs of dates and ranged from 3.4% to 5.4% for the first, depending on aggregation metric, to no significant differences for the second image pair (Table 1). The differences in red to NIR band ratios of the calibration panel reflected pseudo radiance (Tab. SI1) for the three pairs were 0.42, 0.27 and 0.25 (a, b and c in Fig. 5). The mean and median metrics did not show large differences in terms of R2, MAD and bias. R2 values were, however, lower if only the 20% brightest pixels were considered.

| a) Leaf Off | Median | Mean | Mean, 20% brightest |
| R2 | 0.957 | 0.957 | 0.863 |
| MAD | 0.0236 (3.6%) | 0.0226 (3.4%) | 0.0344 (5.4%) |
| Bias (Direct‐Diffuse) | +0.0220 | +0.0211 | +0.0292 |
| t‐test (p‐val) | 0.0047 | 0.0041 | 0.0002 |
| b) Leaf Off | Median | Mean | Mean, 20% brightest |
| R2 | 0.979 | 0.976 | 0.895 |
| MAD | 0.01 (1.5%) | 0.0097 (1.5%) | 0.0203 (3.25%) |
| Bias (Direct‐Diffuse) | −0.0034 | −0.0027 | +0.0017 |
| t‐test (P‐val) | 0.6392 | 0.6974 | 0.8249 |
| c) Leaf On | Median | Mean | Mean, 20% brightest |
| R2 | 0.946 | 0.941 | 0.774 |
| MAD | 0.0096 (1.1%) | 0.0099 (1.2%) | 0.0125 (1.5%) |
| Bias (Direct‐Diffuse) | −0.0086 | −0.0090 | −0.0093 |
| t‐test (P‐val) | 0.0088 | 0.008 | 0.0232 |
For the concrete reference surface, the standard deviation of NDVI over the entire time‐series was 0.0403. A two‐sample t‐test showed that there was no significant difference between the means of the diffuse and direct measurements for this surface (p‐value =0.216, 95 percent confidence interval: [−0.0635, 0.0157]).
NDVI time‐series
Individual crown mean NDVI values of the 20% brightest pixels in the green band were extracted from all drone orthomosaics based on the digitized extents (Fig. 3C). These were combined to represent the mean NDVI and standard deviation time‐series per species within the study site (Fig. 6). There was one observation less for spruce (Picea spp.) due to incomplete coverage of one acquisition, indicated by a gap in the time‐series. Results from an identical analysis which did not include the brightest pixel extraction step and instead aggregated all crown NDVI values are included in the supplementary information (Fig. SI1). The concrete reference surface revealed some significant variability between acquisitions. Some of these effects were visible in crown NDVI time‐series, such as a slight dip and subsequent spike in NDVI at DOYs 157 and 163, most prominent for gorse and goat willow (Salix caprea). Further variability was visible pre green‐up (approx. DOY 80 to 100) which was related to different illumination conditions. This was not as evident when averaging over all crown pixels (Fig. SI1). The reduced contribution of the understorey for the brightest pixels method was expressed as lower NDVI prior to the green‐up (or baseline NDVI) of the overstorey for most species (Fig. 6, Fig. SI1). The mean difference between NDVI pre‐green‐up and post‐green‐up (using means of the first three and last three acquisitions) for deciduous crowns was 0.252 for the brightest pixels method and slightly lower when using the mean of all pixels (0.234).

The time‐series of PAI and mean NDVI for the woodland revealed that the NDVI increased prior to any significant leaf expansion, as PAI remained low (Fig. 7). The NDVI plateaued around DOY 160, whereas PAI continued to rise. This behaviour is reproduced at the individual crown scale (Fig. SI2), though there is more variability in the NDVI‐PAI relationship between species.

Extracted species phenophases
For deciduous species, the logistic function was successfully fitted to most of the individual crown NDVI time‐series (164 out of 168). Fitting failed for four crowns due to noisy NDVI time‐series including outliers, caused by differences in image orthorectification due to complex crown geometry. Mean residual standard errors per species are reported in the supplementary information (Table SI2). The SOS, MOG and SOP DOYs were extracted from the fitted function (Fig. 8). In 14 individual cases (mostly grey willow and sweet chestnut) the SOS could not be extracted due to the first maximum in rate‐of‐change of the curvature of the fitted model lying beyond the bounds of the analysed time‐span (prior to DOY 80) and the corresponding MOG and SOP values were excluded due to poorly constrained fits. The largest within‐species variation in metrics could be observed for SOS, particularly for ash (Fraxinus excelsior), oak and chestnut. A comparison of the same metrics extracted using the mean NDVI of all crown pixels (Fig. SI3) showed that there were differences in green‐up timing depending on the method used. SOS was found to be approx. 1.5 days later, MOG 1 day earlier and SOP 3 days earlier averaged over all classes when using the 20% brightest pixels (Table SI4). Comparisons of the median identified SOS date with mean budburst and first leaf dates of seven species from UKPN observations in South West England revealed that SOS was closer in time to first leaf than to budburst for all species except alder and ash, which were close to budburst (8.4 and 1.11 days difference) (Table SI3). All median SOS dates were within ±1 standard deviation from the mean first leaf dates except for hawthorn and elder where reported dates fall prior to the first drone acquisition.

For evergreen species, the changes in NDVI were very small and the extraction of phenological metrics based on the logistic function are expected to be uncertain. The function could not be fitted to spruce and holly time series due to no consistent NDVI increase. Pines exhibited a small NDVI increase and phenological metrics could be extracted which are reported in the supplementary information (Fig. SI4). Gorse was excluded from this analysis due to flowering‐related variations throughout the spring.
Discussion
Robustness of drone‐based NDVI time‐series
NDVI data from drone‐based multi‐spectral images can effectively track seasonal changes at the crown level. Visual assessments of extracted time‐series separated by species (Fig. 6) showed a clear differentiation between many species with feasible, gradual increase in NDVI for deciduous crowns starting at different times throughout the spring, whereas evergreen crowns such as the spruce showed very little variation.
While optical satellite data has proven effective for classifications of tree species from landscape scale to individual level by incorporating information from different phenological stages or sensors with fine spatial resolution (Immitzer et al., 2019; Persson et al., 2018; Wagner et al., 2018) as well as for monitoring forest phenology (Han et al., 2013; Pastor‐Guzman et al., 2018; Prabakaran et al., 2013; Walker et al., 2012), the latter is predominantly based on medium to coarse resolution data (e.g. Sentinel 2: 10‐20 m, Landsat: 30 m, MODIS: 250 m) which result in mixed pixels for heterogeneous ecosystems such as the study site presented here and are influenced by soil background effects and understorey vegetation (Boyd et al., 2011; Ryu et al., 2014). The ability to resolve individuals and separate species phenology even within small mixed tree stands overcomes these limitations. Drone data, however, include different uncertainties related to platform, acquisition and data processing, which our results suggest should be assessed robustly by users to confirm time‐series reliability.
A prerequisite for the time‐series analysis at the individual crown level is a high georeferencing accuracy of the drone‐based products. Here a horizontal error smaller than one GSD was achieved through the use of 6 well‐distributed GCPs which proved sufficient without requiring further co‐registration of the datasets. The GCP approach remains the method of choice for multi‐temporal studies from drone platforms (Berra et al., 2017; van Iersel et al., 2018), yet in forest environments with canopy occlusion doing so remains a challenge (Wallace et al., 2016). Developments of direct georeferencing using real‐time kinematic (RTK) or post‐processing kinematic (PPK) solutions show promise to facilitate high accuracy drone surveys of these systems (Forlani et al., 2018; Zhang et al., 2019) and marker‐free co‐registration approaches may be applied if absolute geolocation is not necessary to further reduce equipment and survey costs (Poblete et al., 2018; Scheffler et al., 2017).
A further constraint on drone‐based time‐series monitoring is the impact of different illumination conditions on data quality. The effects of dominantly direct or diffuse illumination are challenging to correct (Hakala et al., 2013) and lead to major impacts in band‐ratio indices if individual bands are not perfectly aligned (Berra et al., 2019). The alignment issues were mitigated by the MCA sensor and software‐based workflow utilized in this study as the red and NIR spectral bands were co‐registered during pre‐processing. Assessed over calibration targets, illumination conditions have been shown to have a small impact (differences of 1‐3%) on the reflectance of the red and NIR bands of the Parrot Sequoia (Deng et al., 2018). Differences are expected to be amplified over vegetation due to structural influences on reflectance anisotropy (Sandmeier & Itten, 1999).
We identified a small but significant impact of contrasting illumination on NDVI values for vegetation crowns for a leaf‐off and a leaf‐on direct‐diffuse image pair, whereas a third image pair during leaf‐off conditions showed no significant difference (Fig. 5, Table 1). The different findings for the two leaf‐off image pairs could relate to greater wavelength‐specific differences in irradiance inferred from calibration panel reflected radiance (red to NIR band ratios differed by 0.42 and 0.27 between images) and differences in reconstruction quality of tree crowns. In images of the leaf‐off canopy acquired under direct irradiance conditions, background features were more visible which besides implicitly affecting the radiometry also influenced the DSM and mosaicking algorithm of the photogrammetric software due to issues with crown reconstruction (influencing image projection and prioritizing ground pixels visible within multiple images in some cases). To minimize visible understorey in gaps which can lead to bias in retrieved crown NDVI metrics, acquiring data during overcast conditions is recommended prior to leaf development.
After tree crown green‐up, a small but significant decrease in NDVI could be observed when imaged during direct illumination conditions (bias =ca. −0.01) (Fig. 5, Table 1). This may relate to shading effects within the crown due to its structure, leaf angular distributions and solar illumination angles, which other studies have noted can have a significant effect on index products (Aboutalebi et al., 2019; Damm et al., 2015; Gamon et al., 2006; Schläpfer et al., 2013).
Overall, the identified illumination‐based uncertainties were roughly an order of magnitude smaller (up to 0.034 for leaf‐off, 0.013 for leaf‐on) than the NDVI difference pre‐ and post‐green‐up of deciduous trees (0.252). Therefore, for the study of spring phenology of deciduous crowns, observations from contrasting illumination conditions can be included. Illumination‐based uncertainties must be revisited for the monitoring of photosynthetic phenology in evergreen species as the pigment indices commonly used are much more sensitive to these effects (D’Odorico et al., 2020; Mõttus et al., 2015; Takala & Mõttus, 2016).
Drone‐derived phenology metrics
The high spatial and temporal resolutions that can be achieved with drone‐based NDVI information makes resolving seasonal changes at the individual and species level possible, but their potential for phenological monitoring needs to be assessed critically. Based on the NDVI time‐series, for deciduous crowns, a logistic function could be fitted in most cases and SOS, MOG and SOP metrics derived. Based on a comparison with PAI values for the woodland NDVI was found to mirror the increase in PAI, with some non‐linearity arising from NDVI increase prior to overstorey green‐up due to understorey effects which are known to influence both satellite and in situ sensor‐derived phenology estimates (Ahl et al., 2006; Ryu et al., 2014; Soudani et al., 2012; White et al., 2014), as well as NDVI saturation prior to maximum PAI (Fig. 7). The latter is a previously identified limitation of vegetation indices such as NDVI due to ratioing and the reflectance at red wavelengths experiencing little change over denser canopies (Birky, 2001; Wang et al., 2005; Zhu & Liu, 2015). Changes in the enhanced vegetation index (EVI) have been found to better reflect vegetation phenology (Kowalski et al., 2020) but could not be calculated from Sequoia data due to requiring spectral information at blue wavelengths (Huete et al., 2002). It should be stressed that these conclusions are from satellite‐based studies at coarser resolutions and further work is needed to explore these dynamics in drone‐acquired data. For specific species, the relative timing of the phenological stages matched well with previous observations, such as the early green‐up of alder and hawthorn, the slower progression and late peak greenness of alder and late but rapid green‐up of ash and oak tree crowns (Kuster et al., 2014; Leslie et al., 2017; Mijnsbrugge & Janssens, 2019; Vilhar et al., 2013). On the basis of individual crown PAI measurements (Fig. SI2), differences in green‐up timing between oak and hawthorn could be verified and UKPN observations of species phenology in the same geographic region (counties Cornwall and Devon in South West England) confirmed a late green‐up of beech, ash and oak through mean budburst dates (Table SI3). The median identified SOS date was closer to the date of first leaf emergence and later than budburst for most species, though it should be considered that regional variations of the reported events may occur at smaller spatial scales than the county level. In contrast, Berra et al., (2019) found that SOS identified from drone‐acquired RGB data was generally earlier than SOS derived from field observations (defined as >90% of bud open) in North East England. Differences could arise due to understorey vegetation growth, the vegetation index used or observation timing and fitted model. Here, further observations of pre‐green‐up would enable a more robust model fit and SOS extraction as some species showed increasing NDVI shortly after the first drone acquisitions. The intra‐species variation in phenological metrics were relatively large, particularly for SOS and MOG (Fig. 8), and include outliers. For some species such as oaks, regional differences in photoperiod response have been observed (Basler & Körner, 2012). Erroneous early green‐up and variability may be influenced by the growth of understorey vegetation, particularly bluebells (Hyacinthoides non‐scripta) within the woodland and further uncertainties arose in this experiment from observation gaps in the dataset, caused primarily by adverse weather conditions which did not permit safe drone flights.
From inspecting species‐specific NDVI time‐series (Fig. 6) it was evident that an increase in greenness for evergreen species could also be distinguished. NDVI increased slightly in pine species (Pinus spp.) in early summer which may indicate phenological changes as have previously been observed for pine species using satellite‐based NDVI (Aragones et al., 2019). However, NDVI is not a good indicator of leaf‐level photosynthetic phenology in evergreen species (D’Odorico et al., 2020; Wong et al., 2019).
Flowering events reduced NDVI (by increasing reflectance in both red and NIR) in hawthorn (mid‐May) and evergreen gorse (Bowman et al., 2008). These species‐specific effects could be better represented and analysed with dedicated models and masking, for example the 20% brightest pixel method here showed to be much more sensitive to the hawthorn blooming than the mean of all pixels (Fig. 6, Fig. SI1). Mapping of flowering phenology has also been achieved using drone‐acquired RGB imagery (Neumann et al., 2020) and further comparisons are required to identify which spectral bands can yield the best indicators.
The future of drone‐based phenological studies
In order for the presented methodology to be scalable to larger study areas, some logistical concerns must be addressed. The manual digitizing of individual crowns should ideally be replaced either by object‐based segmentation of crowns and automated species classification making use of the fine spatial resolution textural and multi‐temporal dataset (Xu et al., 2020, example of workflows which could be modified for drone datasets: Guirado et al., 2019; Qiu et al., 2020; Sheeren et al., 2016; Wagner et al., 2018), or auxiliary datasets such as detailed forest inventories, for example based on laser scanning data (Burt et al., 2019; Duncanson et al., 2014). A further limiting factor for the deployment of drones to improve upon phenocam‐based monitoring of phenology is their timely deployment and operational costs considering the fieldwork effort involved for reaching remote sites and the susceptibility to weather conditions. While this study showed that both, direct and diffuse illumination conditions are conducive to the acquisition of useable datasets, gaps in time‐series result from too changeable conditions, too high wind‐speeds and other limitations. For dedicated study sites, fully automated deployment informed, for example, by weather‐station data from a hub including battery recharging station as is being pioneered in robotics and smart‐cities research (Byun et al., 2016; Galkin et al., 2019; Menouar et al., 2017) may soon offer the desired solution.
With sensor technology as well as associated workflows still rapidly evolving, the ideal payload for phenological monitoring has yet to be determined. Saturation and understorey influences as were observed in this study could be addressed by performing leaf‐cover classification at the pixel level (Park et al., 2019). Crucial information on phenology also includes processes such as pigment changes at leaf‐scale, which cannot easily be decoupled from changes in crown leaf area when using merely the NDVI. Including further spectral bands, particularly the red‐edge, in phenological analysis could improve this separation (D’Odorico et al., 2020; Delegido et al., 2011). Finally, new opportunities are arising with the development of drone sensors that exhibit identical bandsets to those of satellite sensors (Revill et al., 2019). These will aid in understanding and validating studies based on coarser‐resolution time‐series from satellites where multiple crowns are aggregated and mixed pixels will show a combination of individual canopy phenology, scaling differently depending on phenological stage (Chen et al., 2018; Klosterman et al., 2018).
Conclusions
Recent studies using drone‐mounted consumer‐grade cameras revealed that tree‐specific phenological monitoring is possible. Here, we showed that the NDVI derived from a calibrated drone‐mounted multi‐spectral sensor can provide the basis for deriving and comparing phenological metrics for many species in a heterogeneous ecosystem. This finely resolved spatial information holds great potential for revealing adaptations in species‐specific phenology to changing local conditions such as the micro‐climate, as well as the possible identification of flowering events with high relevance to pollinator research. The signal from the vegetated understorey, however, exerts an influence, even at the fine spatial resolutions of drone‐based datasets, and should be addressed. Furthermore, we demonstrated that observations from contrasting illumination conditions can introduce bias, but this was relatively small when compared with the NDVI change during spring green‐up. Given conditions throughout each flight are still consistent, data from overcast or clear‐sky acquisitions can be used together for NDVI‐based time‐series analysis which greatly increases the feasibility of these studies in temperate climates with variable weather conditions.
Acknowledgements
Karen Anderson and Dominic Fawcett received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska‐Curie Grant Agreement No 721995. The authors would like to thank Keith Wilson for providing full access to Trelusback farm as a study site.